/*
 * LPSolver.cpp
 *
 * Copyright 2012 Fabian Wenzelmann
 *
 * This file is part of Recognition-Procedures-for-Boolean-Functions.
 *
 * Recognition-Procedures-for-Boolean-Functions is free software: you
 * can redistribute it and/or modify it under the terms of the
 * GNU General Public License as published by the Free Software Foundation,
 * either version 3 of the License, or (at your option) any later version.
 *
 * Recognition-Procedures-for-Boolean-Functions is distributed in the hope
 * that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
 * Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Recognition-Procedures-for-Boolean-Functions.
 *
 * If not, see <http://www.gnu.org/licenses/>.
 */

#include "bfr/lpbLP/LPSolver.hpp"

#include <cassert>

#include <algorithm>
#include <functional>
#include <queue>
#include <utility>
#include <vector>
#include <list>

#include <lp_lib.h>

#include "bfr/Literal.hpp"
#include "bfr/Clause.hpp"
#include "bfr/ClauseSet.hpp"
#include "bfr/WinderMatrix.hpp"
#include "bfr/LPB.hpp"
#include "bfr/lpbLP/DNFTree.hpp"

namespace bfr {

LPSolver::LPSolver(bool setIntModel, tighten_lp t, bool rt)
    : intModel(setIntModel), tighten(t), regularityTest(rt) { }

LPSolver::~LPSolver() { }

LPB* LPSolver::convert(const ClauseSet *cs)
        const throw (LPBConversionException) {
    WinderMatrix wm(cs);
    std::vector<BooleanVector*> mtps(cs->size());
    size_t pos = 0;

    sortMatrix(wm);
    size_t *sortedPos = new size_t[wm.size()];
    // for each variable save the position in the sorted matrix
    // TODO: is there a way to do it smarter?
    for (size_t i = 0; i < wm.size(); i++) {
        sortedPos[wm.getVarIndex(i)] = i;
    }

    // create true points
    // also make a copy of dnf in list representation
    DNFTree::content_t *listRep = new DNFTree::content_t();
    ClauseSet::const_iterator e = cs->end();
    for (ClauseSet::const_iterator it = cs->begin();
            it != e;
            it++) {
        const Clause *c = *it;
        std::list<Literal> *clauseList = new std::list<Literal>;
        BooleanVector *point = new BooleanVector(wm.size(), false);
        mtps[pos++] = point;
        Clause::const_iterator ce = c->end();
        for (Clause::const_iterator cit = c->begin();
                cit != ce;
                cit++) {
            const Literal &l = *cit;
            clauseList->push_back(l);
            const size_t lID = l.getId();
            // should be correct now...
            //point->operator [](lID) = true;
            point->operator [](sortedPos[lID]) = true;
        }
        listRep->push_back(clauseList);
    }

    if (regularityTest) {
        if (!isRegular(listRep, wm, mtps)) {
            for (size_t i = 0; i < mtps.size(); i++) {
                delete mtps[i];
            }
            delete[] sortedPos;
            throw LPBConversionException("DNF is not regular.");
        }
    }

    std::vector<BooleanVector *> mfps;
    computeMFPs(mtps, mfps);

    lprec *lp = setupLP(mtps, mfps, wm, intModel, tighten);
    //print_lp(lp);
    try {
        LPB *solution = solveLP(lp, wm);
        // delete the minimal true points
        for (size_t i = 0; i < mtps.size(); i++) {
            delete mtps[i];
        }
        for (size_t i = 0; i < mfps.size(); i++) {
            delete mfps[i];
        }
        delete_lp(lp);
        delete[] sortedPos;
        return solution;
    }
    catch (const LPBConversionException &e) {
        for (size_t i = 0; i < mtps.size(); i++) {
            delete mtps[i];
        }
        for (size_t i = 0; i < mfps.size(); i++) {
            delete mfps[i];
        }
        delete[] sortedPos;
        delete_lp(lp);
        throw e;
    }
}

void LPSolver::sortMatrix(WinderMatrix &wm) {
    size_t variables = wm.size();
    std::sort (wm.getMatrix(),
               wm.getMatrix() + variables,
               CompareMatrixEntryGe(variables));
}

bool LPSolver::isRegular(const DNFTree::content_t *dnf, const WinderMatrix &matrix,
                         const std::vector<BooleanVector* > &mtps)
    const {
    if (matrix.size() == 0)
        return true;

    // create the DNFTree
    DNFTree tree;
    // create the root node entry
    // we know that it is not false because the matrix has size > 0
    tree.createRoot(dnf, false);

    // create a queue that stores the node ids that must be explored
    std::queue<size_t> waiting;
    // push the id of the first node, this is 0
    waiting.push(0);
    {
        std::pair<const DNFTree::content_t*, const DNFTree::content_t*>
            splitResult(0, 0);
        bool firstFinal;
        bool secondFinal;
        size_t depth;
        size_t varIndex;
        // process all entries
        while (!waiting.empty()) {
            size_t next = waiting.front();

            if (tree.getFinal(next)) {
                // we don't want to split it
                // just delete the entry
                assert (next != 0);
                tree.deleteEntry(next);
                waiting.pop();
                continue;
            }

            // so now it is not final so we have to split
            depth = tree.getDepth(next);
            assert (depth < matrix.size());
            const DNFTree::content_t *nextDNF = tree.getDNFEntry(next);
            varIndex = matrix.getVarIndex(depth);
            splitAway(nextDNF, varIndex, splitResult,
                    &firstFinal, &secondFinal);

            // TODO TODO:
            // this should be possible to do
            // because each node is only processed once
            // but it seems that one node is popped more than
            // once?
            // AND the delete only destroys the DNF, but the error
            // occurs when setting the right child?
            //tree.deleteEntry(next);

            const DNFTree::content_t *first = splitResult.first;
            const DNFTree::content_t *second = splitResult.second;

            if (firstFinal) {
                if (first->size() != 0) {
                    size_t leftID =
                            tree.createLeftChild(next, first, true);
                    waiting.push(leftID);
                }
            }
            else {
                size_t leftID =
                        tree.createLeftChild(next, first, false);
                waiting.push(leftID);
            }

            if (secondFinal) {
                if (second->size() != 0) {
                    size_t rightID =
                            tree.createRightChild(next, second, true);
                    waiting.push(rightID);
                }
            }
            else {
                size_t rightID =
                        tree.createRightChild(next, second, false);
                waiting.push(rightID);
            }
            waiting.pop();
        }
    }

/*    for (size_t i = 0; i < tree.size(); i++) {
        std::cout << "Node with DNF " << *tree.getDNFEntry(i) << std::endl
                  << "Children are Nodes " << tree.getLeftChild(i) << ", "
                  << tree.getRightChild(i) << std::endl;
    }*/


    // after creating the tree we start the regularity test
    const size_t numRuns = matrix.size() - 1;

    for (size_t i = 0; i < numRuns; i++) {
        for (size_t k = 0; k < mtps.size(); k++) {
            BooleanVector *mtp = mtps[k];
            if ((!mtp->operator [](i)) && (mtp->operator [](i + 1))) {
                // change the positions in the point, after the implicant test
                // we will change them again
                mtp->operator [](i) = true;
                mtp->operator [](i + 1) = false;
                if (!isImplicant(tree, mtp, matrix)) {
                    return false;
                }
                mtp->operator [](i) = false;
                mtp->operator [](i + 1) = true;
            }
        }
    }
    return true;
}

void LPSolver::splitAway(const DNFTree::content_t *dnf, size_t variableID,
        std::pair<const DNFTree::content_t*,
                  const DNFTree::content_t*> &result,
        bool *firstFinal, bool *secondFinal) {
    *firstFinal = false;
    *secondFinal = false;
    // TODO maybe estimate size somehow?
    DNFTree::content_t *first = new DNFTree::content_t();
    DNFTree::content_t *second = new DNFTree::content_t();
    DNFTree::content_t::const_iterator e = dnf->end();
    //const Literal currentLiteral(variableID);
    for (DNFTree::content_t::const_iterator it = dnf->begin();
         it != e;
         it++) {
        const std::list<Literal> *c = *it;
        bool containsVariable = false;
        std::list<Literal> *newClause = new std::list<Literal>();
        for (std::list<Literal>::const_iterator cit = c->begin();
             cit != c->end();
             cit++) {
            const Literal &l = *cit;
            if (l.getId() == variableID) {
                // we set the flag that we found the variable
                // but we don't add it
                assert (!containsVariable);
                containsVariable = true;
            }
            else {
                // add it to the new clause
                newClause->push_back(l);
            }
        }
        if (containsVariable) {
            first->push_back(newClause);
            if (newClause->size() == 0)
                *firstFinal = true;
        }
        else {
            second->push_back(newClause);
            if (newClause->size() == 0)
                *secondFinal = true;
        }
    }
    if (first->size() == 0) {
        *firstFinal = true;
    }

    if (second->size() == 0) {
        *secondFinal = true;
    }
    result.first = first;
    result.second = second;
}

bool LPSolver::isImplicant(const DNFTree &tree,
                           const BooleanVector *mtp,
                           const WinderMatrix &matrix) const {
    size_t u = 0;
    for (size_t k = 0; k < mtp->size(); k++) {
        if (tree.isLeaf(u))
            return true;

        size_t leftChild = tree.getLeftChild(u);
        size_t rightChild = tree.getRightChild(u);

        if (mtp->operator [](k)) {
            if (leftChild != DNFTree::npos) {
                u = leftChild;
                continue;
            }
            else {
                assert (rightChild != DNFTree::npos);
                u = rightChild;
            }
        }
        else {
            if (rightChild != DNFTree::npos) {
                u = rightChild;
                continue;
            }
            else {
                return false;
            }
        }
    }
    assert (tree.getLeftChild(u) == DNFTree::npos);
    assert (tree.getRightChild(u) == DNFTree::npos);
    return true;
}

void LPSolver::computeMFPs(std::vector<BooleanVector *> &mtps,
            std::vector<BooleanVector *> &result) {
    std::sort(mtps.begin(), mtps.end(), LPSolver::ComparePoints());
    size_t *nu = new size_t[mtps.size()];
    // TODO necessary?
    for (size_t i = 0; i < mtps.size(); i++) {
        nu[i] = 0;
    }
    size_t vars;
    bool val1, val2;
    for (size_t i = 1; i < mtps.size(); i++) {
        vars = mtps[i]->size();
        for (size_t j = 0; j < vars; j++) {
            val1 = mtps[i - 1]->operator [](j);
            val2 = mtps[i]->operator [](j);
            if ((!val1) && val2) {
                nu[i] = j + 1;
                break;
            }
        }
    }

    for (size_t i = 0; i < mtps.size(); i++) {
        BooleanVector *point = mtps[i];
        vars = point->size();
        for (size_t j = nu[i]; j < vars; j++) {
            //if (mtps[i]->operator [](j)) {
            if (point->operator [](j)) {
                assert (nu[i] <= j);
                BooleanVector *newPoint = new BooleanVector(*point);
                newPoint->operator [](j) = false;
                for (size_t k = j + 1; k < vars; k++) {
                    newPoint->operator [](k) = true;
                }
                result.push_back(newPoint);
            }
        }
    }

    delete[] nu;
}

lprec* LPSolver::setupLP(const std::vector<BooleanVector *> &mtps,
                        const std::vector<BooleanVector *> &mfps,
                        const WinderMatrix &wm,
                        bool intModel,
                        tighten_lp tighten) throw (LPBConversionException) {
    const size_t numVariables = wm.size();
    const size_t cols = numVariables + 1;
    lprec *result = make_lp(0, cols);

    if (result == 0) {
        throw LPBConversionException("Unable to create LP.");
    }

    set_add_rowmode(result, true);

    REAL *row = new REAL[cols];
    int *colno = new int[cols];

    // first of all add non-negative constraint for all variables
    // if intModel is true also make it a int variable
    // we don't have to add the non-negative constraint for the variables
    // because lpsolve sets the lower bounds to 0 by default
    if (intModel) {
        for (size_t i = 0; i < numVariables; i++) {
            row[0] = 1;
            colno[0] = i + 1;
            if (!set_int(result, i + 1, true)) {
                delete[] row;
                delete[] colno;
                delete_lp(result);
                throw LPBConversionException("Unable to mark variable as"
                        " integer (out of memory?).");
            }
        }
        if (!set_int(result, cols, true)) {
            delete[] row;
            delete[] colno;
            delete_lp(result);
            throw LPBConversionException("Unable to mark variable as"
                    " integer (out of memory?).");
        }
    }
    /*for (size_t i = 0; i < numVariables; i++) {
        row[0] = 1;
        colno[0] = i + 1;
        if (intModel) {
            // mark variable as int
            if (!set_int(result, i + 1, true)) {
                delete[] row;
                delete[] colno;
                delete_lp(result);
                throw LPBConversionException("Unable to mark variable as"
                        " integer (out of memory?).");
            }
        }
        if (!add_constraintex(result, 1, row, colno, GE, 0)) {
            delete[] row;
            delete[] colno;
            delete_lp(result);
            throw LPBConversionException("Unable to add a constraint"
                    " to the LP (out of memory?).");
        }
    }

    // TODO is this correct? should be...
    if (intModel) {
        if (!set_int(result, cols, true)) {
            delete[] row;
            delete[] colno;
            delete_lp(result);
            throw LPBConversionException("Unable to mark variable as"
                    " integer (out of memory?).");
        }
    }*/

    for (size_t i = 0; i < mtps.size(); i++) {
        BooleanVector *point = mtps[i];
        size_t pos = 0;
        // create a new column
        for (size_t j = 0; j < numVariables; j++)  {
            if (point->operator [](j)) {
                // we add a constraint for this variable
                row[pos] = 1;
                colno[pos] = j + 1;
                pos++;
            }
        }
        row[pos] = -1;
        colno[pos] = cols;
        pos++;
        if(!add_constraintex(result, pos, row, colno, GE, 0)) {
            delete[] row;
            delete[] colno;
            delete_lp(result);
            throw LPBConversionException("Unable to add a constraint"
                                " to the LP (out of memory?).");
        }
    }

    for (size_t i = 0; i < mfps.size(); i++) {
        BooleanVector *point = mfps[i];
        size_t pos = 0;
        // create a new column
        for (size_t j = 0; j < numVariables; j++) {
            if (point->operator [](j)) {
                row[pos] = 1;
                colno[pos] = j + 1;
                pos++;
            }
        }
        row[pos] = -1;
        colno[pos] = cols;
        pos++;
        // TODO
        // lp solve does not have a < constraing, only <=
        // so is this correct?
        if (!add_constraintex(result, pos, row, colno, LE, -1)) {
            delete[] row;
            delete[] colno;
            delete_lp(result);
            throw LPBConversionException("Unable to add a constraint"
                                " to the LP (out of memory?).");
        }
    }

    // check if we should add some extra constraints and also check what type
    // of constraints
    if (tighten == tighten_lp::NEIGHBOURS) {
        // only add a constraint for neighbouring variables, that
        // means for variable v_i and variable v_{i+1}
        // so we have to compare all neighbouring variables
        // we know that v_i >= v_{i+1}
        // so we only have to check if they are equal
        LPSolver::CompareMatrixEntry comp(numVariables);
        size_t **matrix = wm.getMatrix();
        // create an array of size 2 that we can use to add the constraint
        // we always set one factor to one and the other factor to
        // -1 and add the constraint that it must be equal to 0
        REAL currRow[2] = {1, -1};

        int currColNo[] = {-1, -1};
        for (size_t i = 1; i < numVariables; i++) {
            // check if variable an position i - 1 is equal to
            // variable on position i
            bool added = false;
            currColNo[0] = i;
            currColNo[1] = i + 1;
            short compRes = comp(matrix[i - 1], matrix[i]);
            assert (compRes >= 0);
            if (compRes == 0) {
                // they are equal so add an equals constraint
                added = add_constraintex(result, 2, currRow, currColNo, EQ, 0);
            }
            else {
                // just add an >= constraint
                added = add_constraintex(result, 2, currRow, currColNo, GE, 0);
            }
            if (!added) {
                delete[] row;
                delete[] colno;
                delete_lp(result);
                throw LPBConversionException("Unable to add a constraint"
                                    " to the LP (out of memory?).");
            }
        }
    }

    else if ((tighten == tighten_lp::ALL) && (numVariables > 0)) {
        // we want to create a constraint for all variable pairs
        // idea is the same as above
        LPSolver::CompareMatrixEntry comp(numVariables);
        size_t **matrix = wm.getMatrix();
        REAL currRow[] = {1, -1};
        int currColNo[] = {-1, -1};
        // TODO: are the loops correct?
        // little bit tricky
        std::vector<bool> eqResults;
        eqResults.reserve(numVariables - 1);

        // in the first iteration, when we create all entries
        // for v_0, we also initialize the values of eqResults
        size_t i = 1;
        // we compare the matrix entry by entry, this variable is used
        // to signal that the first variable is still equal to the current
        // entry
        bool stillEqual = true;
        for (size_t j = i; j < numVariables; j++) {
            currColNo[0] = i;
            currColNo[1] = j + 1;
            short compResult = comp(matrix[j - 1], matrix[j]);

            int constrType = -1;
            assert (compResult >= 0);

            if (compResult == 0) {
                eqResults.push_back(true);
                if (stillEqual) {
                    constrType = EQ;
                }
                else {
                    // we found equal variables but no longer equal to
                    // v_{i - 1}
                    constrType = GE;
                }
            }
            else {
                eqResults.push_back(false);
                stillEqual = false;
                constrType = GE;
            }
            assert ((constrType == GE) || (constrType == EQ));
            // now add constraint
            if (!add_constraintex(result, 2, currRow, currColNo, constrType, 0)) {
                delete[] row;
                delete[] colno;
                delete_lp(result);
                throw LPBConversionException("Unable to add a constraint"
                                    " to the LP (out of memory?).");
            }
        }
        // now add all other constraints
        // its just a look up in the results vector
        i = 2;
        while (i < numVariables) {
            size_t j = i;
            currColNo[0] = i;
            // first find the longest sequence s. t. the variables are equal
            while((j < numVariables) && (eqResults[j - 1])) {
                currColNo[1] = j + 1;
                if (!add_constraintex(result, 2, currRow, currColNo, EQ, 0)) {
                    delete[] row;
                    delete[] colno;
                    delete_lp(result);
                    throw LPBConversionException("Unable to add a constraint"
                                        " to the LP (out of memory?).");
                }
                j++;
            }
            // now add all >= constraints
            while (j < numVariables) {
                currColNo[1] = j + 1;
                if (!add_constraintex(result, 2, currRow, currColNo, GE, 0)) {
                    delete[] row;
                    delete[] colno;
                    delete_lp(result);
                    throw LPBConversionException("Unable to add a constraint"
                                        " to the LP (out of memory?).");
                }
                j++;
            }
            i++;
        }
    }

    delete[] row;
    delete[] colno;

    return result;
}

LPB* LPSolver::solveLP(lprec *lp, const WinderMatrix &matrix)
        throw (LPBConversionException) {
    size_t numVariables = matrix.size();
    set_add_rowmode(lp, false);
    set_verbose(lp, CRITICAL);
    // TODO work with return value
    std::cout << "Before solve(lp)" << std::endl;
    int convResult = solve(lp);
    std::cout << "After solve(lp)" << std::endl;
    if (convResult != FEASFOUND && convResult != OPTIMAL) {
        throw LPBConversionException("There is no feasible "
                "solution for the LP.");
    }
    size_t cols = numVariables + 1;
    REAL *result = new REAL[cols];
    if (!get_variables(lp, result)) {
        delete[] result;
        throw LPBConversionException("Error getting the solution via"
                " get_variables.");
    }
    //for (size_t i = 0; i < cols; i++) {
    //    std::cout << result[i] << std::endl;
    //}
    LPB *lpb = new LPB(result[cols - 1], numVariables);
    for (size_t i = 0; i < numVariables; i++) {
        lpb->addSummand(result[i], Literal(matrix.getVarIndex(i)));
    }
    delete[] result;
    return lpb;
}

} // namespace bfr
